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»(t)  - F(t)8(t)  + P(t)H"(t)(y(t)  - H(t)S(t))  , 

■ ' ' x(0)  - 0 (2a 

where  P(t)  , the  covariance  of  the  error  of  the 
state  estimate,  obeys  the  Riccati  equation 

P(t)  - F(t)P(t)  + P(t)F''^(t)  + G(t)C^(t)  - 

- P(t)H’'(t)H(t)P(t) 

P(0)  - n (2b 


Abstract 


A simple  differential  equation  for  the  triangular 
square-root  of  the  error  covariance  of  the  linear 
state  estimator  is  derived.  Previous  algorithms 
involved  an  antisymmetric  matrix  in  the  square- 
root  differential  equation.  In  the  constant 
model  case,  Chandrasekhar- type  equations  are 
shown  to  constitute  a set  of  fast  square-root 
algorithms  for  the  derivative  of  the  error  vari- 
ance. 

Square-root  algorithms  for  the  smoothing 
problem  are  presented  and  as  in  the  discrete  case, 
an  array  method  for  handling  continuous  square- 
roots  is  developed.  This  method  also  yields 
very  naturally  the  usual  normalizations  of 
stochastic  calculus,  suggesting  extensions  to 
more  general  stochastic  equations,  even  to  esti- 
mators for  nonlinear  models. 


Over  the  years,  the  Riccati  equation  has  been 
studied  in  great  detail,  and  as  one  way  to  improve 
its  nunerical  conditioning,  a whole  family  of 
alternative  square-root  algorithms  has  been 
introduced.  In  the  discrete-case,  square-root 
algorithms  have  been  studied  by  Potter  [2], 

Golub  [3],  Schmidt  [4],  Kaminski  and  Bryson  [5], 
(6],.Blerman  [7],  [8],  Morf  and  Kallath  (9], 
among  others. 

Continuous-time  square-root  algorithms  have 
attracted  somewhat  less  interest  (see  Andrews 
[10],  Tapley  and  Choe  [11],  Blerean  [12],  and 
Morf  and  Kallath  [9],  Appendix  C).  Somewhat 
surprisingly,  all  known  solutions  explicitly 
introduce  a certain  antisymmetric  matrix  into  the 
differential  equation  for  the  square-roots.  In 
Section  II  of  this  paper,  we  present  a differen- 
tial equation  for  the  triangular  square-roots 
that  does  not  explicitly  contain  any  such  anti- 
symnetrlc  matrix.  We  show  that  these  matrices 
are  generators  of  the  .orthogonal  transformations 
that  relate  the  various  square-roots  and  that  a 
particular  choice  of  this  matrix  will  result  in 
triangular  square-roots.  Among  other  results,  we 
shall  also  discuss  the  stability  of  our  new 
equations,  and  we  shall  also  present  some  alterna- 
tive forms  (such  as  the  information  filter  forms 
for  the  case  of  high  initial  uncertainty  of  the 
state  estimate).  In  passing  we  also  note  that 
these  ideas  can  also  be  applied  to  the  Chan- 
drasekhar-type algorithms  (for  constant  models). 


I.  Introduction 


Since  a full  version  of  this  paper  will  appear 
in  [1],  we  present  here  only  an  outline  of  our 
results. 

We  assume  that  we  are  given  a state-space 
model 

x(t)  - F(t)x(t)  + G(t)u(t) 

y(t)  - H(t)x(t)  + v(t)  (la) 

where  u(>)  and  v(*)  are  white  noises  with 
normalized  Joint  covariance: 


Here,  T denotes  the  transpose  and  E the  expec- 
tation. ^ 

Then,  x(t) , the  linear  least  squares  estimate 
of  x(t)  given  the  observed  data 

" {y(8)  , 0 < s < t)  can  be  obtained  via  the 
o — — 

Kalman  filter  equation 


II.  Continuous-time  square-root  algorithms 


If  S is  a square-root  of  P and  P(*)  la 
strictly  positive  definite,  so  that  we  can  write 
P • SS^  where  S is  nonsingular,  then  the 
Riccati  equation  (2b)  can  be  rewritten  as 

P-SS^  + SsT-(F-IiS  sTh'^h)s  + 

+ »1  C G^  s“^s'^  + S sT(fT  - Lh^H  S S^)  + 

+ ^ S S-1  G . (3) 

This  equation  is  satisfied  if  the  square-root 
S obeys  the  differential  equation 
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(10) 


S - (F  - >1  S H)S  + >S  G S“^  (Aa) 

with  Initial  conditions 

S(0)  - . <*b) 

However,  It  la  well-known  that  matrix  square- 
roots  are  not  unique,  for  If  S Is  a square-root, 
so  la  ST  , where  T Is  any  orthonormal  matrix. 
To  study  how  this  non-uniqueness  Is  reflected 
Into  the  differential  equation,  we  note  first 
that  any  differentiable  orthonormal  matrix  T(t} 
can  be  uniquely  characterized  by  a skew-or-anti- 
synmetric  matrix  generator  a(-)  such  that 

d(*)+a^(*)-0  (5a) 


T(t)  - A(t)T(t) 


T - I 
o o n 


Therefore,  if  S is  the  square-root  obeying 
equation  (4)  and  If 

- ST  (6a) 

(L 

by  differentiation  of  (6a),  we  obtain  the  dif- 
ferential equation  satlsfl^  by  S^ 

S - (F  - **S  sVh)S  (S  a's^  + »1GG^)S"'^  , 
o o o ^ A . 

S-(0)  -ST.  '(6b) 

<l  o o 

In  this  equation,  the  aiatrlx 
A(t)  A S(t)a(t)sT(t)  is  antlsymetrlc  and  can 
be  made  arbitrary  since  C(t)  Is  arbitrary. 
Therefore,  this  establishes  that  all  the  square- 
roots  of  P are  given  by  the  solutions  of  the 
differential  equation 

S « (F  - JJSS'^H^H)S  + (A  + »jGg'')S~’' 


S(0)  - S, 


s s^  - n 

o o o 


idiere  A(*)  Is  an  arbitrary  antlsynnetrlc 
matrix. 

Clearly,  the  square-root  matrix  S solution 
of  (7)  depends  of  the  antlsynnetrlc  matrix 
A(-)  chosen.  This  dependence  gives  us  an  addi- 
tional degree  of  freedom,  which  can  be  exploited 
in  different  ways.  For  example,  a difficulty 
%rf.th  formula  (7)  Is  that  S~l  has  to  be  computed 
at  each  integration  step.  However,  Andrews  (10] 
noted  that  A(*)  could  be  chosen  so  as  to  make 
S lower  triangular,  so  that  S~^  can  be  com- 
puted store  easily.  To  do  this,  rewrite  (7)  as 

S - (A  + r)s“'^  («) 

with  r F SS^  + ^(GG*  - SsVhSS^). 

Then,  S can  be  maintained  In  lower  triangular 
form  by  calculating  A such  that 

t(A  + r)S~^l_  - 0 

and  by  Inserting  this  value  of  A In  (8).  (Here 
[Ml_  denotes  the  strictly  upper  triangular  part 
of  the  matrix  N ) . 

A somewhat  more  direct  approach  was  discovered 
by  Tapley  and  Choe  [11]  who  noted  that  (8),  when 
rewritten  as 


A + r - S S (10) 

was  giving  a lowcr-t lines  uppcr-trlangular  ' 
factorisation  of  A -f  T ,,  and  ther^y  were  able  to 
simultaneously  compute  S (In  lower  triangular 
form)  and  A . 

Actually,  It  Is  possible  to  find  a differential 
equation  for  S In  lower  triangular  form  that 
does  not  explicitly  Involve  any  antisymmetric 
matrix  In  the  equation:  multiply  (3)  on  the  left 
by  S”^  and  on  the  right  by  S"~  to  obtain 

S”^S  + - M 

where  M M + f'*'  + cc''^  - h’'^H 

and  F - S"*FS  , G - s"^G  , H - HS  . 

-1* 

Since  S Is  lower  triangular,  S S la  the 
lower  triangular  part  of  M and  we  can  identify 
(an  idea  analogous  to  the  Wlener-Hopf  technique) 

S - SlMl^/2  • 

Here,  ( the 'lower-triangular  part" 

operator  defined  by 


^"U/2ij  ■ i’*“ll 


for  1 < j 
.for  1 - j 
for  1 < j 


for  an  arbitrary  matrix  U . The 'Strictly  lower- 
triangular  part' has  zero  values  for  1 - J . 
Similarly,  ( • the  "upper-triangular  part" 

operator  is  defined  by 

W1-/2  • . 

The  Slain  advantage  of  equation  (11)  over  those 
of  Andrews  and  Tapley  and  Choe  Is  chat  It  does 
not  Involve  explicitly  any  antlsymnetrlc  matrix 
In  the  differential  equation.  Moreover,  the 
matrices  (F  , .G,  H)  have  a nice  Interpretation; 
these  Biatrlces  arise  in  Che  dynamic  model  for  Che 
modified  state-variable  n - S~^x  . The  signifi- 
cance of  this  choice  'of  variables  is  that  the 
variance  of  the  error  n will  be  1 . The 
dynamic  model  for  n can  be  obtained  from  (1) 
and  (11)  as 

n - (F  - (M]^^2)n  + Cu  , y-Hn  + v (12) 

Stability  of  the  equation 

Square-root  methods  have  several  advantages 
over  the  Rlccatl  equation  (2b):  for  example, 

P - SS^  Is  always  non-tiegatlve  definite,  a 
property  which  is  not  always  guaranteed  by  direct 
numerical  Integration  of  (26).  Furthermore,  since 
the  "condition  nua^er"  of  S is  the  square-root 
of  the  one  for  F , F can  be  computed  with 
greater  accuracy.  On  the  other  hand,  these 
algorithms  require  a slightly  larger  asiounc  of 
computations,  as  Is  shown  In  Appendix  B of  [1] 
where  computational  aspects  are  discussed. 

It  should  also  be  pointed  out  that  the  sta- 
bility properties  of  Che  Rlccatl  equation  are 
conserved  by  equation  (11):  as  In  [14],  assume 
that  the  pair  (F,G)  (respectively  (H,F))  Is 
uniformly  completely  controllable  (respectively. 


observable)  and  that  Sj^  and  $2  are  two  solu-* 
Lions  of  the  square-root  equation  (11)  correspond- 
Ing'to  different  Initial  conditions  and 

S20.,.duch  that  SjoSJq  “ Hjq  > 0 and 

®5o®io  " ^20  ^ 0 • _ f 

In  this  case,  Pj  ■ S^SJ  and  P2“  S2S2  •'f« 
positive-definite  solutions  of  the  Riccatl 
equation  (2b)  and  if  dP  ■ P^  - P?  , it  was 
proved  by  Kalman  in  (14]  Chat  dP(t)  •*  0 as 

t -*  • . 

Now,  denote  AS  " Sj  - $2  , then 
AP  • AS  . S^  -t-  $2  . As'^  and  by  sniltiplication  on 
the  left  by  Sjl  and  on  the  right  by  Sj'^  , we 
obtain 

S“^AS  + AS’^S"^  - S"^AP  Sj^  . (13) 

But  $2  and  83  are  both  lower  triangular, 
and  therefore  ST^/^S  is  the  lower  cciangular  part 
of  SJ^APJ^: 

AS  - SjlS~^AP  S2'''^]^^2 

and  since  AP(t)  -►  0 ■ as  t “ and  S2  and 
are  bounded  for  i “ 1,2  (P2  and  are 

proved  to  be  bounded  in  [14]  we  can  conclude 
that  AS(t)  0 as  t . 

This  Indicates  that  Che  square-root  equation 
(11)  and  the  Riccatl  equation  (2b)  have  the  same 
stability  properties. 

Information  filter  forma 


In  several  Instances,  for  example  when 
(the  uncertainty  on  the  initial  condition  x(0)) 
becomes  large,  it  is  more  convenient  to  compute 
P“l(t)  instead  of  P(t)  . Filters  involving 
P“1  or  its  square-root  JJ(P~1  ■ 0 have 
been  called  information  filters. 

In  this  case,  f)  obeys  the  differential 
equation 

n . (-f"^  - GG^)n  + (a  + »sH^H)n“'^  (is) 


where  n(0)tl'^(0)  • and 
antisymmetric  matrix. 

Then,  Instead  of  computing  x 
convenient  to  compute  3 ■ P“^  x 
information  filter 


is  an  arbitrary 

, it  is  more 
with  the 


3 - -(F  + CG^Ol^)^  d + H^y  . (16) 

T 

Alternatively,  we  can  also  propagate  n - 0 
via  the  variance  normalized  information  filter 

8-(a->*G5'^->j  H^n  + H‘'^y  (17) 

where  o ^ » G ■ fi'^C  , H « , 

We  have  seen  in  (12)  that  such  a filter  arises 
naturally  in  the  Interpretation  of  equation  (11) 
for  the  triangular  square-root  S . The  informa- 
tion filter  form  of  equation  (11)  is 

iP  - “IMJ+/2 


can  be  obtained  by  propagating  the  square-root  of 
f(t)  Instead  of  the  square-root  of  P(t)  . How- 
ever our  earlier  discussion  suggests  that  we 
should  consider  a family  of  square  roots,  differ- 
ing from  each  other  by  orthogonal  transformations. 

To  see  how  this  can  be  Incorporated  let  us 
review  Che  derivation  of  Che  Chandrasekhar  equa- 
tions. The  first  step  is  to  exploit  the  constancy 
of  the  model  parameters  by  differentiating  the 
Riccatl  equation  (2b)  to  obtain 

P - (F  - K(t)H)P  + P(F  - K(t)H)^  (19) 

and  therefore  P(t)  - ^Vc(ttO)P(0)i>J(t,0) 

where  ♦i,(t,0)  is  the  transition  matrix  associ- 
ated to  (F  - K(t)H)  , K(t)  - P(t)HT  being  the 
Kalman  gain. 

Now,  consider  for  simplicity  the  special  case 
where  llg  >,0  (known  initial  conditions).  In 
(his  case  P(0)-  GGI  , so  that  we  can  factor 
P(t)  ■ L(t)LT(t)  in  square-root  form,  t(t) 
being  nxo  with  o - rank  GC^  £ q:  the  number 
of  inputs.  . 

Then  substituting  P “ LL^  in  (19),  we  find 
that  L(c}  obeys  the  Chandrasekhar-type  equations 
(13] 

Ut)  - (F  - K(t)H)L(t)  - L(t)a(t) 

K(t)  - L(t)LT(t)HT  (20) 

where  L(0)lT(0)  - GGT,.  K(0)  - n<,HT  , a(t)  being 
some  axa  antisymmetric  matrix. 

In  (12],  a(t)  ••  0 was  chosen,  however  the 
introduction  of  a(t)  in  (20)  can  provide  some 
additional  control  over  the  numerical  behavior 
of  the  square-root  L(t)  . 

III.  Conclusions 

In  [1],  we  show  that  the  results  of  Section  II 
can  be  extended  to  Che  fixed-point  smoothing 
problem.  Also,  continuous-time  array  methods 
similar  to  Chose  developed  in  [9]  by  Horf  and 
Kailath  for  the  discrete-time  case  are  presented. 
This  approach  yields. very  naturally  the  usual 
normalizations  of  stochastic  calculus,  suggesting 
extensions  to  more  general  stochastic  equations. 

In  Appendix  A of  [1],  we  discuss  an  alternative 
set  of  variance  normalized  Chandrasekhar-type 
equations  that  involve  only  one  differential 
equation  for  Che  square-root  of  the  error 
covariance;  however,  the  variance  normalized  state 
model  parameters  have  to  be  obtained.  In 
Appendix  B,  we  present  operation  counts  for  Che 
algorithm  of  Section  II  and  those  of  Andrews  and 
Tapley  and  Choe.  The  square-root  algorithms 
involve  roughly  15  to  302  more  computations  Chan 
the  Riccatl  equation.  Our  new  equation  has  the 
advantage  of  giving  an  explicit  differential 
equation  for  Che  square-root,  which  is  consequently 
easier  to  implement  for  analog  simulation. 


where  0 is  the  upper  triangular  square-root  of 
• 

Time-Invariant  Systems. 

In  the  constant  model  case  (H,  F,  G constanOi 
it  is  known  [13]  that  fast  square-root  algorithms 
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